# Get earliest VA DC 
# Enrolled for at least a year and not dead prior
# 2.3M
db.handle <- odbcDriverConnect('server=vhacdwa01.vha.med.va.gov; database=CDWWork; driver={SQL server}; trusted_connection=true')

cohort  <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #cohort
SELECT
  SSN_NBR, RATING_RCVD_DT
INTO #cohort
FROM(
  SELECT 
    SSN_NBR, RATING_RCVD_DT, ROW_NUMBER() OVER(PARTITION BY SSN_NBR ORDER BY RATING_RCVD_DT ASC) AS row_num
  FROM OMHSP_VBA.App.RITA_2021_23663 AS vba
  INNER JOIN CDWWork.SPatient.SPatient AS pat
    ON vba.SSN_NBR=pat.PatientSSN AND (DeathDateTime>RATING_RCVD_DT OR DeathDateTime IS NULL)
  INNER JOIN CDWWork.Patient.Enrollment AS enroll
    ON pat.PatientSID=enroll.PatientSID AND enroll.EnrollmentDate<=DATEADD(day, -480, RATING_RCVD_DT)
  WHERE RATING_RCVD_DT>=CAST('2004-01-01' AS date) AND RATING_RCVD_DT<CAST('2018-01-01' AS date) AND EnrollmentDate<CAST('2016-07-30' AS date)
) AS z
WHERE row_num=1

SELECT * FROM #cohort")

cohort <- data.table(cohort)
cohort$Year <- as.numeric(substring(cohort$RATING_RCVD_DT, 1, 4))
cohort$RATING_RCVD_DT <- as.Date(cohort$RATING_RCVD_DT)

inpat <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;
SELECT
  SSN_NBR, RATING_RCVD_DT, RelativeMonth, COUNT(DISTINCT VisitDate) AS N_Inpat, Inpat=1
FROM(
SELECT SSN_NBR, RATING_RCVD_DT, CAST(AdmitDateTime AS date) AS VisitDate, FLOOR(DATEDIFF(day, RATING_RCVD_DT, AdmitDateTime)*1.0/30) AS RelativeMonth
FROM CDWWork.Inpat.Inpatient AS inpat
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON inpat.PatientSID=pat.PatientSID
INNER JOIN #cohort
  ON pat.PatientSSN=#cohort.SSN_NBR AND ABS(DATEDIFF(day, RATING_RCVD_DT, AdmitDateTime))<=510 AND NOT RATING_RCVD_DT=CAST(AdmitDateTime AS DATE)
WHERE AdmitDateTime>=CAST('2002-07-01 00:00:00' AS datetime2(0)) AND AdmitDateTime<CAST('2020-01-01 00:00:00' AS datetime2(0))
) AS z
GROUP BY SSN_NBR, RATING_RCVD_DT, RelativeMonth
")
inpat$RATING_RCVD_DT <- as.Date(inpat$RATING_RCVD_DT)

outpat <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; DROP TABLE IF EXISTS #stopcode;
SELECT DISTINCT StopCodeSID, CASE WHEN StopCode=130 THEN 1 ELSE 0 END AS ED, CASE WHEN StopCode=450 THEN 1 ELSE 0 END AS CP INTO #stopcode FROM CDWWork.Dim.StopCode;

SELECT 
  SSN_NBR, RATING_RCVD_DT, FLOOR(DATEDIFF(day, RATING_RCVD_DT, VisitDateTime)*1.0/30) AS RelativeMonth, 
  COUNT(DISTINCT CAST(VisitDateTime AS date)) AS N_Outpat,
  MAX(CASE WHEN (psc.ED=1 OR ssc.ED=1) THEN 1 ELSE 0 END) AS ED
FROM CDWWork.Outpat.Workload AS visit
LEFT OUTER JOIN #stopcode AS psc
  ON visit.PrimaryStopCodeSID=psc.StopCodeSID
LEFT OUTER JOIN #stopcode AS ssc
  ON visit.SecondaryStopCodeSID=ssc.StopCodeSID 
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON visit.PatientSID=pat.PatientSID
INNER JOIN #cohort
  ON pat.PatientSSN=#cohort.SSN_NBR AND ABS(DATEDIFF(day, RATING_RCVD_DT, VisitDateTime))<=510 AND NOT RATING_RCVD_DT=CAST(VisitDateTime AS DATE)
WHERE VisitDateTime>=CAST('2002-07-01 00:00:00' AS datetime2(0)) AND VisitDateTime<CAST('2020-01-01 00:00:00' AS datetime2(0))
AND (psc.CP=0 AND ssc.CP=0)
GROUP BY SSN_NBR, RATING_RCVD_DT, FLOOR(DATEDIFF(day, RATING_RCVD_DT, VisitDateTime)*1.0/30)
")
outpat$RATING_RCVD_DT <- as.Date(outpat$RATING_RCVD_DT)


opioids <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off;
SELECT 
  SSN_NBR, RATING_RCVD_DT, FLOOR(DATEDIFF(day, RATING_RCVD_DT, ReleaseDateTime)*1.0/30) AS RelativeMonth, Opioid=1
FROM CDWWork.RxOut.RxOutpatFill AS rx
INNER JOIN OMHSP_PERC_Share.DOEx.Lookup_MorphineEquivalence AS perc
  ON rx.NationalDrugSID=perc.NationalDrugSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON rx.PatientSID=pat.PatientSID
INNER JOIN #cohort
  ON pat.PatientSSN=#cohort.SSN_NBR AND ABS(DATEDIFF(day, RATING_RCVD_DT, ReleaseDateTime))<=510 AND NOT RATING_RCVD_DT=CAST(ReleaseDateTime AS DATE)
WHERE Opioid IS NOT NULL
AND ReleaseDateTime>=CAST('2002-07-01 00:00:00' AS datetime2(0)) AND ReleaseDateTime<CAST('2020-01-01 00:00:00' AS datetime2(0))
GROUP BY SSN_NBR, RATING_RCVD_DT, FLOOR(DATEDIFF(day, RATING_RCVD_DT, ReleaseDateTime)*1.0/30)
")
opioids$RATING_RCVD_DT <- as.Date(opioids$RATING_RCVD_DT)


nsaids <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off
DROP TABLE IF EXISTS #nsaids;

SELECT DISTINCT NationalDrugSID, DrugNameWithoutDose INTO #nsaids FROM CDWWork.Dim.NationalDrug AS a
INNER JOIN CDWWork.Dim.DrugClass AS dim
  ON a.PrimaryDrugClassSID=dim.DrugClassSID
INNER JOIN CDWWork.Dim.DrugNameWithoutDose AS b
  ON a.DrugNameWithoutDoseSID=b.DrugNameWithoutDoseSID
WHERE DrugClassCode IN ('CN104', 'MS102') AND NOT DrugNameWithoutDose='*Unknown at this time*';

SELECT 
  SSN_NBR, RATING_RCVD_DT, FLOOR(DATEDIFF(day, RATING_RCVD_DT, ReleaseDateTime)*1.0/30) AS RelativeMonth, NSAID=1
FROM CDWWork.RxOut.RxOutpatFill AS rx
INNER JOIN #nsaids
  ON rx.NationalDrugSID=#nsaids.NationalDrugSID
INNER JOIN CDWWork.SPatient.SPatient AS pat
  ON rx.PatientSID=pat.PatientSID
INNER JOIN #cohort
  ON pat.PatientSSN=#cohort.SSN_NBR AND ABS(DATEDIFF(day, RATING_RCVD_DT, ReleaseDateTime))<=510 AND NOT RATING_RCVD_DT=CAST(ReleaseDateTime AS DATE)
WHERE ReleaseDateTime>=CAST('2002-07-01 00:00:00' AS datetime2(0)) AND ReleaseDateTime<CAST('2020-01-01 00:00:00' AS datetime2(0))
GROUP BY SSN_NBR, RATING_RCVD_DT, FLOOR(DATEDIFF(day, RATING_RCVD_DT, ReleaseDateTime)*1.0/30)
")
nsaids$RATING_RCVD_DT <- as.Date(nsaids$RATING_RCVD_DT)

mortality <- sqlQuery(db.handle, "set nocount on; set ansi_warnings off; 
SELECT #cohort.SSN_NBR, CAST(MIN(DeathDateTime) AS Date) AS DeathDate FROM #cohort INNER JOIN CDWWork.SPatient.SPatient AS pat ON #cohort.SSN_NBR=pat.PatientSSN GROUP BY #cohort.SSN_NBR
")

mortality <- merge(x=mortality, y=cohort[,c("SSN_NBR", "RATING_RCVD_DT")], by.x="SSN_NBR", by.y="SSN_NBR")
mortality$DeathDate <- as.Date(mortality$DeathDate, "%Y-%m-%d")
mortality$RATING_RCVD_DT <- as.Date(mortality$RATING_RCVD_DT, "%Y-%m-%d")
mortality$DeathRelativeMonth <- floor(as.numeric(mortality$DeathDate-mortality$RATING_RCVD_DT)/30)

# Create panel dataset
panel <- data.table(SSN_NBR=rep(cohort$SSN_NBR, each=33), RelativeMonth=rep(-16:16, length(unique(cohort$SSN_NBR))));
panel <- merge(x=panel, y=cohort, by.x="SSN_NBR", by.y="SSN_NBR")
panel <- merge(x=panel, y=inpat, by.x=c("SSN_NBR", "RATING_RCVD_DT", "RelativeMonth"), by.y=c("SSN_NBR", "RATING_RCVD_DT", "RelativeMonth"), all.x=T)
panel <- merge(x=panel, y=outpat, by.x=c("SSN_NBR", "RATING_RCVD_DT", "RelativeMonth"), by.y=c("SSN_NBR", "RATING_RCVD_DT", "RelativeMonth"), all.x=T)
panel <- merge(x=panel, y=opioids, by.x=c("SSN_NBR", "RATING_RCVD_DT", "RelativeMonth"), by.y=c("SSN_NBR", "RATING_RCVD_DT", "RelativeMonth"), all.x=T)
panel <- merge(x=panel, y=nsaids, by.x=c("SSN_NBR", "RATING_RCVD_DT", "RelativeMonth"), by.y=c("SSN_NBR", "RATING_RCVD_DT", "RelativeMonth"), all.x=T)
panel <- merge(panel, y=mortality[,c("SSN_NBR", "RATING_RCVD_DT", "DeathRelativeMonth")], by.x=c("SSN_NBR", "RATING_RCVD_DT"), by.y=c("SSN_NBR", "RATING_RCVD_DT"), all.x=T)
panel$Dead <- (panel$RelativeMonth>=panel$DeathRelativeMonth); panel$Dead[which(is.na(panel$Dead))] <- 0
setnafill(panel, cols=c("N_Inpat", "Inpat", "N_Outpat", "ED", "NSAID", "Opioid", "Dead"), fill=0)
panel$Outpat <- (panel$N_Outpat>0)
#rm(inpat, nsaids, opioids, outpat, pain)

panel$Cohort <- format(panel$RATING_RCVD_DT, "%Y-%m")
panel$CalendarMonth <- format(panel$RATING_RCVD_DT+panel$RelativeMonth*30, "%Y-%m")
panel$CalendarYear <- substring(panel$CalendarMonth, 1, 4)
panel$RelativeMonth <- as.factor(panel$RelativeMonth)
panel$RelativeMonth <- relevel(panel$RelativeMonth, "-12")

panel <- panel[(DeathRelativeMonth>16 | is.na(DeathRelativeMonth))]

memory.limit(size=100000)
N_outpat <- summary(feols(N_Outpat~RelativeMonth|CalendarYear+Cohort, data=panel), cluster="Cohort")$coeftable
Outpat <- summary(feols(Outpat~RelativeMonth|CalendarYear+Cohort, data=panel), cluster="Cohort")$coeftable
ED <- summary(feols(ED~RelativeMonth|CalendarYear+Cohort, data=panel), cluster="Cohort")$coeftable
Inpat <- summary(feols(Inpat~RelativeMonth|CalendarYear+Cohort, data=panel), cluster="Cohort")$coeftable
Opioid <- summary(feols(Opioid~RelativeMonth|CalendarYear+Cohort, data=panel), cluster="Cohort")$coeftable
NSAID <- summary(feols(NSAID~RelativeMonth|CalendarYear+Cohort, data=panel), cluster="Cohort")$coeftable


N_outpat <- N_outpat[,c("Estimate", "Std. Error")]; N_outpat <- rbind(N_outpat[1:4,], c(0,NA), N_outpat[5:32,]); N_outpat <- data.table(N_outpat); N_outpat$Estimate <- N_outpat$Estimate + mean(panel$N_Outpat[which(panel$RelativeMonth=="-12")])
Outpat <- Outpat[,c("Estimate", "Std. Error")]; Outpat <- rbind(Outpat[1:4,], c(0,NA), Outpat[5:32,]); Outpat <- data.table(Outpat); Outpat$Estimate <- Outpat$Estimate + mean(panel$Outpat[which(panel$RelativeMonth=="-12")])
ED <- ED[,c("Estimate", "Std. Error")]; ED <- rbind(ED[1:4,], c(0,NA), ED[5:32,]); ED <- data.table(ED); ED$Estimate <- ED$Estimate + mean(panel$ED[which(panel$RelativeMonth=="-12")])
Inpat <- Inpat[,c("Estimate", "Std. Error")]; Inpat <- rbind(Inpat[1:4,], c(0,NA), Inpat[5:32,]); Inpat <- data.table(Inpat); Inpat$Estimate <- Inpat$Estimate + mean(panel$Inpat[which(panel$RelativeMonth=="-12")])
Opioid <- Opioid[,c("Estimate", "Std. Error")]; Opioid <- rbind(Opioid[1:4,], c(0,NA), Opioid[5:32,]); Opioid <- data.table(Opioid); Opioid$Estimate <- Opioid$Estimate + mean(panel$Opioid[which(panel$RelativeMonth=="-12")])
NSAID <- NSAID[,c("Estimate", "Std. Error")]; NSAID <- rbind(NSAID[1:4,], c(0,NA), NSAID[5:32,]); NSAID <- data.table(NSAID); NSAID$Estimate <- NSAID$Estimate + mean(panel$NSAID[which(panel$RelativeMonth=="-12")])

output <- cbind(N_outpat, Outpat, ED, Inpat, Opioid, NSAID)
names(output) <- c("outpatCount_coef", "outpatCount_se", "outpat_coef", "outpat_se", "ed_coef", "ed_se", "inpat_coef", "inpat_se", "opioid_coef", "opioid_se", "nsaid_coef", "nsaid_se")
rownames(output) <- c("EventMonth-16", "EventMonth-15", "EventMonth-14", "EventMonth-13", "EventMonth-12", "EventMonth-11", "EventMonth-10", "EventMonth-9", "EventMonth-8", "EventMonth-7", "EventMonth-6", "EventMonth-5", "EventMonth-4", "EventMonth-3", "EventMonth-2", "EventMonth-1",
                      "EventMonth0", "EventMonth1", "EventMonth2", "EventMonth3", "EventMonth4", "EventMonth5", "EventMonth6", "EventMonth7", "EventMonth8", "EventMonth9", "EventMonth10", "EventMonth11", "EventMonth12", "EventMonth13", "EventMonth14", "EventMonth15", "EventMonth16")
fwrite(output, "O:/OMHSP_PERCAnalyt/Disability Ratings/Disability Event Study/va_eventstudy_application.csv")


